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TECHNICAL MEMORANDUM 


ON THE DESIGN OF STRUCTURAL COMPONENTS USING MATERIALS 
WITH TIME-DEPENDENT PROPERTIES 

I. INTRODUCTION 


With the increased use of polymer matrix composite materials for aerospace applications, 
design engineers are faced with the need for knowledge of environmental effects on the functional 
behavior of structures. A variety of epoxies are widely used as a binding agent or matrix for fiber 
reinforced composites. Thermoset polymers are often used over thermoplastic polymers because of 
their better thermal stability and chemical resistance. A great advantage of thermoset polymers is 
their higher resistance to creep and stress relaxation. Typical thermoset matrix materials are 
epoxies, polyesters, and vinyl esters. Although they have relatively better creep characteristics than 
thermoplastics, they will still “relax” as a function of time, load, temperature, and other environ- 
mental factors. The understanding of this viscoelastic behavior when designing with composite 
materials is the topic of this report. 

During the typical preliminary design phase of composite material structures, basic classical 
lamination theory (CLT) solutions are used to obtain stress, strain, deflections, natural frequency, 
and buckling strength of the structure analyzed. With the increased emphasis on long-term 
aerospace structures (10 to 30 years useful life), it is of great importance that viscoelastic effects 
also be included in the early design phase in order to obtain knowledge of the structural/functional 
behavior of the hardware after an extended period of time. With analytical representation of the 
material characteristics, the composite design could possibly require modification of initial dimen- 
sions and geometry in order to meet critical functional requirements at the end of their useful life. 

This report presents one of the widely used methods, namely the elastic-viscoelastic corre- 
spondence principle, in the analysis of viscoelastic structural materials. For this report, the concern 
is with the basic time dependency of the viscoelastic material properties. The goal is to present the 
design engineer with an effective method to analyze, in closed form or numerically, and to design 
time-dependent structures. The problem investigated is one of the most basic problems in structural 
design, “beams.” 

It is the author’s goal to provide information for structural design engineers to consider during 
the inception and preliminary design of any structure with time-dependent properties. The methods 
presented are proven and widely accepted yet simple to understand and apply. 


II. THE ELASTIC BEAM PROBLEM 


The classical problem of pure bending of beams with constant cross section is chosen, since 
the basic equations for stress, strain, deflection, and natural frequency have been well established. 1-3 
For purposes of this report, a cantilever beam with a single load applied at its free end will be con- 
sidered. The solution derivation will have the usual assumptions. These are: 


1. The beam is thin. This implies that the thickness is much smaller than any of the other 
physical dimensions. 


2. The deflection of the beam in the direction of the applied load is small compared to the 
beam thickness. (This assumption has been shown to be applicable even in the case of 
relatively large resulting deflections.) 4 

3. The in-plane strains e x , By, and e ^ are small compared to unity. 

4. Transverse shear strains e yi and e xl are negligible. 

5. The material obeys Hooke’s law. 

6. Rotatory inertia terms are negligible. 

7. There are no body forces. 

8. The material is isotropic. 


A. Determination of Stress and Strain 

^ The bending stress in a one-dimensional beam can be determined from the following equa- 


Mc 

°b=— ( 1 ) 

where M is the bending moment, c is the distance from the neutral axis to the outermost surface, 
and / is the moment of inertia of the beam cross section. For the beam configuration shown in 
figure 1, the bending moment can be expressed as: 

M = -P(L-x), (l a ) 

where 


P = pb 


c - -t/2 (for the top (tension) surface of the beam) 

1 = bt 3 /ll (for a rectangular cross section). 
Substituting these values into equation (1), one obtains: 

6P(L-x) 


( 2 ) 


2 


. .IN 




Figure 1. Cantilever beam configuration. 

Following Hooke’s, law one can express the axial strain caused by the applied bending 
moment by simply dividing equation (2) by the modulus of elasticity (Young s modulus) of the beam 
material. In this manner one obtains: 

6 P(L-x) rT i 


B. Determination of Deflection 

The formulation which links the curvature of the central line of the beam cross section with 
the applied bending moment is called the Euler-Bernoulli law. This is expressed in differential form 
as: 

d 2 w M /, 


With the knowledge that the slope and deflection at the fixed end of the beam are zero, equation (4) 
can be integrated twice to yield the following expression for the deflection at the free end of the 
beam: 


C. Determination of Natural Frequency 

For the one-dimensional beam problem, classical linear elastic beam theory yields the 
following frequency equation in the absence of in-plane forces: 


_ r d 4 w d 2 w 

El y + m T7 

dx 4 dt 2 


= 0 . 
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With the knowledge that the slope and deflection vanish at the fixed end and the shear force and 
moment are zero at the free end, this equation can be used to obtain the various bending vibration 
modes of the beam. The resulting characteristic equation of the beam is: 


where 


cos XL cosh XL -- 1 , 
2 

^4 _ mca 

El ’ 


(7) 

( 8 ) 


Solving equation (7) for X L and substituting the results into equation ( 8 ) yields the following 
expression for the first or natural frequency of the beam: 


co = 3.51601 


El 
ml L 4 


(9) 


m. THE SINGLE-PHASE VISCOELASTIC BEAM PROBLEM 


A. Material Properties Characterization 

In order to obtain a solution to any structural design problem where the material analyzed 
has time-dependent properties such as plastics, elastomers, and resin-based matrix composite 
materials, an analytical expression for the relaxation modulus E(t) and Poisson’s ratio v(t) must be 
developed. The relaxation modulus is obtained from standard stress relaxation tests. In these tests, 
the viscoelastic material is subjected to a constant strain. Under the influence of this strain, the 
material will relax; and the stress will gradually decrease. The stress is measured at specific time 
intervals, and the relaxation modulus is plotted as a function of time. These data are then curve fitted 
to a function which can be readily manipulated to perform the necessary analysis for the solution of 
the problem. A very common function used is the Prony series curve fit due to Gaspard Francois 
Clair Marie Riche de Prony (1755-1839). A form of this function is: 


fit) = A + £s,e r ' f . ( 10 ) 

i=i 

Recent investigations have produced methods for obtaining the coefficients and exponents of this 
function automatically with the aid of computers . 5 6 For purposes of this report, the following function 
will be used for the relaxation modulus: 


E{t) = A+B l e-rUB 2 e~r2‘+B 3 e-^ t , ( 1 1 ) 

where 

A = 180,000 
B l = 5,000 
B 2 = 5,000 
S 3 = 170,000 


4 


ft = 1,000 

y 2 = 10 

ft = 0.10 • 


The time-dependent Poisson’s ratio is obtained from the relaxation modulus data and from 
the knowledge of the bulk modulus of elasticity of the material. The bulk modulus of elasticity can be 
expressed as: 


m 

3(1 - 2 v(0) ' 


( 12 ) 


Solving for v(t ) one obtains: 



m 

6K 


(13) 


Substituting equation (11) into equation (13) and expanding, one obtains the following expression: 

v(t) = A n +B n e-^'+B 23 e-^ t +B 33 e-ri' , (14) 

where 


A _ 1 A 

A i 3 - — dj 3 - 


B \ 3 ~ “ 


2 6 K 

B i 


6K 


B 33 “ “■ 


*2. 

6K 

i± 

6K 


Figures 2 and 3 show plots of the relaxation modulus and Poisson’s ratio for the Prony functions of 
equations (11) and (14), respectively. 



Log of time (hours) 

Figure 2. Relaxation modulus of a single-phase viscoelastic material. 
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Figure 3. Poisson’s ratio of a single-phase viscoelastic material. 


B. Elastic-Viscoelastic Correspondence Principle 

A very effective method of solution of time-dependent structural problems is obtained by 
applying the Laplace transform to the time-dependent functions. This removes the time variable, and 
the analysis problem for the viscoelastic body is converted to an “associated elastic” problem. This 
method allows the solution of the viscoelastic problem by simply expressing the constitutive equa- 
tions and boundary conditions as functions of the Laplace transform parameter s . Once a solution to 
the “associated elastic” problem is obtained in the Laplace domain, it can be inverted into the origi- 
nal time domain, and a solution to the original viscoelastic problem is developed. For a single-phase 
material exhibiting properties that are time dependent, the solution of the beam problem becomes 
fairly straightforward. In fact, by using the “associated elastic” problem approach which is officially 
known as the “elastic-viscoelastic correspondence principle,” this problem and many others can be 
solved in closed form. This method has been successfully used and documented by many authors. 7-9 
It will be used in this rep ort to demonstrate its simplicity and usefulness to the design engineer 
when designing structures with time-dependent materials. 

The solutions derived in this report assume linearly viscoelastic materials. This means that, 
for any time interval, the time-dependent functions can be assumed proportionally linear to the 
applied constant load. Also a material is assumed linearly viscoelastic if the combined effects of two 
or more simultaneously applied loads or displacements can be expressed as the sum of the individ- 
ual effects when the same loads or displacements are applied separately. 
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Equation (11) can be expressed in the Laplace domain as: 

— A B, Bj B-i /ic\ 

E(s)= - + L + — 2- + — . ( 15 > 

5 5 + 7j S+Y 2 s+Y 3 


The “associated elastic” expression for the relaxation modulus can be obtained by multiplying 
equation (15) by the Laplace parameter sJ 8 In this manner, one obtains: 


E(s) = sE(s) = 


A , ^1 ! sB 2 , sB 2 . 

s + Yi s+Y 2 s+ Y3 


(16) 


Following this same procedure, the expression for the associated elastic Poisson s ratio is 
expressed as: 


v(s) = sv(s) = 


A + i%_ + i£23_ + ^33_ 
5 + 7! 5+7 2 5+7 3 


(16a) 


C. Determination of Time-Dependent Strain 

Any applied loads, either constant or varying with time, must also be transformed into the 
Laplace domain. So if, in the time domain, equation (3) is expressed as: 


where 


e b (t) = 


sixmo 

E(t) 


_ 6(L-x) (where t is thickness) , 
bt 2 


(17) 


then, in the Laplace domain, equation (17) becomes: 


e b is) = 




Pq 

2 E(s)_ ' 


(18) 


Substituting equation (16) into (18) and combining all terms, one obtains, after some algebraic 
manipulations: 


where 


£ b (s) = S(x)P 0 


\ 5 3 +/i5 2 +/ 2 5 + / 3 ^ 

+ L.25^ + L 3 5 2 + L^S j 


(19) 


7 


fi = 7i + 72 + 73 

h= 7i + 73 

— A + + B 2 + B$ 

f 2 ~ 7\7i + 7 2 73 + 737\ 

fi = 7 x 73 

L 2 = Af x + B\f 4 + B 2 f 6 + Sj/g 

£ 

£ 

11 

< 

> 

11 

+ 

L 3 ~ A f 2 + B xfs + # 2/7 + ^ 3/9 

f A = 72 + 73 

f 9 = 7x72 

II 

> 


fs= 72 73 

The ratio of polynomials in equation (19) can be expressed as: 

■ y3 + /i y2 + hj + h _ F a( s ) ^20) 

s(l x s 3 + L 2 s 2 + L 3 s + L 4 ) sG a (s ) 

The roots of the cubic equation G fl (s) can be obtained, allowing the expression to be written as: 

G fl (5) = (5 + A 1 )(5 + A 2 )(j + A 3 ) . (21) 

It is important to point out that for the physical problem, these roots should never have an imaginary 
component. This knowledge can be used as a check to verify that the numerical calculations have 
been performed appropriately. In fact, all roots in equation (21) should be real and negative for a 
material with a modulus that can be characterized as exponentially decaying. 

One should notice that equation (20) is a quotient of two polynomials with no common fac- 
tors, and the degree of the numerator is lower than that of the denominator. This is the classical 
fraction that can be solved in a straightforward manner by application of Heaviside’s partial fraction 
expansion. 10 Following Heaviside’s procedure, the total derivative of the denominator of equation 
(20) can be expressed as: 


or 


fK«] = s i2M + M GaU) , 

as ds ds 


— |sG fl (s)l = + 3L 2 s 2 + 2L 3 s+ La . 

ds " ' 


Equation (20) can now be expressed as a sum of partial fractions as: 


where 


W _ M x , M 2 f M 3 ( M 4 
sG a Cs) s 5 + Aj s + X 2 s + X 3 


Mi 


w 


S=0 


( 22 ) 

(23) 


(24) 
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JJIllliliil 


" 2 = -i 


Mj = 


W 


-[ SG „( S )] 




M 3~ U 


w 


-[ SC „( S )1 


-A, 


4 ~ A 

ds 


W 


[*<*(*>] 


• = A, 


Substituting equation (24) into equation (19) yields the expression for axial strain due to bending in 
the Laplace domain: 


£ b (s) = S(x)P 0 


M, M 2 M 3 M a 

— + , + , + , 

S 5 + Aj S + A 2 S + A 3 j 


(25) 


Taking the inverse transformation of equation (25) yields the final expression for the time-depen- 
dent axial strain due to bending: 


E b (t) = S(x)P 0 i [ M l +M 2 e X '‘ + M 3 e kl ' + M 4 e . 


(26) 


D. Determination of Time-Dependent Deflection 

Equation (5) is the expression for maximum deflection of the cantilever beam under study. 
Following the same procedure as for the viscoelastic strain calculations, the deflection can be writ- 
ten as a function of the time-dependent variables: 


MO = Q 


p{0 

E(t) 


(27) 


where 


4 if 

Q = — *■ (where t is thickness) 

br 


In the Laplace domain, equation (27) can be written as: 


w(s)= Q\ 


s 2 E(s) 


(28) 
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One should notice that the only difference between equation (28) and equation (18) is in the 
replacement of the term S(x ) with the constant Q. Since both these terms are independent of time, it 
is a straightforward matter to express the time-dependent deflection as in equation (25). This 
expression is: 


w(t) = QP 0 ^M l + M 2 e X, ‘ + M 3 e Xl '+M 4 e . 


(29) 


IV. THE TWO-PHASE VISCOELASTIC BEAM PROBLEM 


The formulation for strain and deflection obtained in section III is applicable for linearly 
viscoelastic materials which are isotropic. Although many plastics can be analyzed in this manner, 
the majority of the viscoelastic materials used for aerospace structural applications are two phase. 
This means that they are composed of one phase that exhibits time-dependent properties and 
another phase that does not. This is the case for many composite materials, in particular the epoxy- 
or resin-based two-phase systems. For example, graphite/epoxy, boron/epoxy, and silicon- 
carbide/epoxy are considered two-phase composites. 


A. Micromechanics Determination of Material Properties 

In order to identify the time-dependent material properties of the two-phase composite, we 
must look at the interaction between the time-dependent and time-independent components at a 
microscopic level. This heterogeneous look at the composite system is known as micromechanics. In 
this report, the classical stiffness approach to micromechanics is used. 11 There are some basic 
restrictions on the composite material. For example, the composite ply (lamina) resulting from the 
constituent parts must be macroscopically homogeneous and macroscopically orthotropic. It must be 
linearly viscoelastic and initially stress-free. For the constituents, the fibers are homogeneous, 
linearly elastic, isotropic, regularly spaced, and perfectly aligned; the matrix is homogeneous, linearly 
viscoelastic, and isotropic. In addition, the bonds between the fibers and the matrix are assumed to 
be perfect (no voids). Although these restrictions are seemingly stringent, modern manufacturing 
methods combined with material characterization at a macroscopic level (£ n , E 2 2> G 12 , etc ) can 

used to “back-out” the necessary constituent characteristics {E m (t),Ej, v m (r), efc.). 

The binder or matrix of the composite material used in this report has a relaxation modulus 
described by equation (11). In this manner, one has: 

E m (t) = A l + B l e- r ' t + B 2 e- r ' , + B 3 e- Y '' . (30) 

Following the micromechanics approach to stiffness, a “rule of mixtures” expression for the appar- 
ent time-dependent Young’s modulus in the direction of the fibers can be obtained. This is: 

E n (t)=E f V f +E m (,)V m , (31) 
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where 


Ef = Young’s modulus for an isotropic fiber 
V j - fiber volume content for the composite material 


E m (t) = relaxation modulus for the matrix 


V m = matrix volume content for the composite material. 

In the direction transverse to the fibers, the apparent Young’s modulus is expressed as: 


£22(0 = 


E f V m +E m (t)V f ' 


(32) 


Several approaches for the accurate determination of the apparent in-plane shear modulus G 12 U) 
have been investigated. Using a variational analysis approach, Foye 12 developed an expression for a 
square array of fibers in the laminate. This represented the best closed-form estimates of this 
orthotropic constant for a unidirectional composite ply. It is expressed, in this report, as follows: 


G i 2 (f) 


GJt) (4-g)+7rAT(r) ( 4N(t) 

2 4 (4-n)N(t)+7c\ ’ 


where 


G m (t) = 


EJt ) 


2[l+v„W] 


W(<) = 


G f (x+4V f ) + G m (t)(x-4Vf) 
G f (x-4V f ) + G m U){x + 4V r ) 



Poisson’s ratio for the matrix material will be described by equation (14) or: 

v m (t) = A 1 3 + B l3 e~ r ' t + B^e-^+B^e-^ . 


(33) 

(34) 

(35) 

(36) 

(37) 


The major Poisson’s ratio for the unidirectional composite lamina can be written using the rule of 
mixtures in the same manner as the time-dependent Young’s modulus E n (t): 


v n (t)=v f V f + v m (t)V m , 


(38) 
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where 


vy = Poisson’s ratio for isotropic fibers. 

The minor Poisson’s ratio is defined as usual from the symmetric properties of the compliance 
matrix: 

• (39) 

EllM 

In order to obtain the time-dependent strain values of the individual plies within the laminated beam, 
one can use the elastic viscoelastic correspondence principle. Due to the fact that the determination 
of strains in composite laminates is critically dependent on stiffness parameters for each ply, as well 
as stiffness parameters for the laminate, a closed-form solution to the time-dependent strains is 
considerably more elaborate than what is expressed in equation (26). In fact, one will realize that 
even for the simplest problems (cantilever laminated beam), although a closed-form solution is 
possible, it is not time nor cost effective to expect a design engineer to obtain them. For more 
complex mathematical models, the function to be inverted is often known only for discrete positive 
real values of the transform parameter therefore making it very difficult if not impossible to obtain an 
exact solution. A more effective approach is the application of numerical inversion methods to obtain 
the approximate transformed solution. A widely used and effective method of inversion is the 
collocation technique due to Schapery. 13 


B. The Collocation Method 

This numerical inversion technique is readily applicable to a general class of problems that 
have a solution of the form: 9 

f(t)=r l +r 2 t+g(t) , (40) 

where Fj and r 2 are constants, and g(t) is the transient component of the solution. The transient 
component is normally expressed, approximately, as a sum of exponential functions or: 

m 

£(0 = ZV' /a ’ > (41) 

V— 1 

where h v and a v are constants. 

The time-dependent axial strain due to bending can be written according to equations (40) 
and (41) as: 

m 

e b (t) = T l +r 2 t + J i h v e~ ,/a - . (42) 

V=1 

After the material experiences the creep that is characteristic of viscoelastic materials, it is assumed 
that the long-term value of strain is approximately constant. In reality, this long-term strain is not 
constant, but for many materials, its rate of change is very small. With this assumption, the linear 
time-dependent component of equation (42) vanishes, yielding: 


12 


imp i 1 



( 43 ) 


e»W = r, + f ]h v e-' /a - . 

V=1 

In the Laplace domain, equation (43) can be expressed as: 

se b (s) = sT l (s) + sg(s) , 


(44) 


where 


and 




1 ’ 


I'M = I 

V=1 


K 

(S + 1/CC V ) ' 


Equation (44) is now written as: 


se b (s) = Tj 


+ 


f s Jh> 

h (^v«v) ' 


(45) 


In order to try to minimize the error of the approximation given by g{s), the transform of the approx- 
imate solution should be equal to the transform of the exact solution, at least at the m discrete 
values of s : 


approx 


\s=ya, ’ 


si?) exact I = 8 (*) 

where 

w = l,2,3,...m . 

Considering equation (46), equation (45) can be expressed as: 

m sh 

se b (s) = T l + ^j- 7 — -r (w = 1,2,3..., m) 

V=1 (V®w"^V®v) 

At t = 0, equation (43) can be solved for the constant T v It is written as: 


(46) 


(47) 


r 1 - 


V=1 


(48) 
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Letting s = l/ocw on the right-hand side of equation (47) and substituting equation (48) into equation 
(47) yields, after rearranging: 


m 


K 


"(1 + aja w ) 


= e b (t 0 )-se b (s) 


\s=ya. 


The expression for a j is: 9 

(Xj = e {J = 1,2,3 . 
Substituting equation (50) into the left-hand side of equation (49) yields: 


I 


K 

l + e 2(w-v)- 


= £b{t 0 )-se b ( s ) 


\s=ya. 


(49) 


(50) 


(51) 


All quantities in the set of equations (51) are either known or can be determined at the discrete 
values of the Laplace parameter 5 except the constants h v . Solving for the h v 's gives the necessary 
information to evaluate equation (48). A final substitution into equation (43) yields the expression 
for the time-dependent strain. 


C. Determination of Time-Dependent Strain 

For the problem of a multilayered viscoelastic composite beam, the time-dependent response 
is readily obtained using an appropriate numerical method. The collocation method is used in this 
section to obtain the axial strain of the cantilever beam due to bending. 


For symmetric laminates under bending loads, the stresses in the 
expressed as: 


<*x 

(*) 

011 012 Q\6 

(*) 

Yx 


. = 

Qn Ql 2 026 

< 

r y ' 

, a *y. 


.016 026 066. 


Yxy 


k ,h ply of the beam can be 


(52) 


where Q { j are the transformed reduced stiffnesses of the ply and y are the curvatures. 11 14 From the 
moment-curvature constitutive relations for bending of composite laminates, one can write: 


£ 

^ 


"a*i a* 2 aV 


M x 

Yy 

, - 

* * * 
A 2 D 22 D 26 

< 

M y 

y*y. 


A*6 »26 066. 


M *y. 
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where D,* are the components of the inverted bending stiffness matrix. For one-dimensional beam 
problems, the following assumption is made: 

M y = Mjy = 0 . (54) 


Expressions for the transformed reduced stiffnesses can be found in Jones 11 and Whitney 14 and will 
not be repeated here. Substituting equations (53) and (54) into (52) yields the expressions relating 
stresses to the ply and laminate stiffnesses, and the applied moment. For the axial direction 
(maximum bending stress direction), the stress is expressed as: 





D n +Qn D\i + Qvb 



From the classical lamination theory, the bending stiffnesses are expressed as: 


J jt=i 


4 



(55) 


(56) 


Once the stresses have been determined, one can express the strains in terms of the stresses by 
transformation of the strain-stress relations from principal material directions to body coordinates. 
The resultant expression is: 



(*) 







(*) 

r 1 

e x 

5ii 

Su 

5l6 

Ox 

£y 

. = 

5 12 

S 2 2 

$26 


Oy > 

. £ *y. 


/l6 

$26 

^66 _ 


°xy 


(*) 


(57) 


where S y are the components of the transformed compliance matrix for the k th ply. Once again 
expressions for the components of the transformed compliance matrix can be found in Jones 11 and 
Whitney 14 and will not be repeated here. 

Using the elastic-viscoelastic correspondence principle, the solution for the time-dependent 
strains can be obtained. The procedure is as follows: 

1. Determine the analytical expressions for the time-dependent relaxation modulus and 
Poisson’s ratio (equations (30) and (37)). 

2. Obtain the Laplace transforms of these functions and determine the associated elastic 
expressions (equations (16) and (16a)). 

3. Calculate the values of the longitudinal and transverse properties for the unidirectional ply 
by transforming equations (31), (32), (33), (38) and (39) into the Laplace domain and 
substituting the values from step 2. 
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4. Calculate the reduced stiffness matrix, compliance matrix, transformed reduced stiffness 
matrix, and transformed compliance matrix in terms of the unidirectional ply properties 
from step 3. 

5. With the information in step 4, calculate the laminate bending stiffness matrix (equation 
(56)). 

6. Identify the applied loads (moments) in the Laplace domain. For a constant moment, this 
is simply dividing the moment by the Laplace parameter. 

7. Calculate stresses and strains using equations (52) and (57). 

8. Express the calculated strains as the “associated elastic” solution by simply multiplying 
the calculated strains from step 7 by the Laplace parameter s. 

9. Solve for the constants, h v , from equation (51). In expanded form, equation (51) can be 
written as: 

1 1 

2 1 + <T 2 

1 1 


l + e 2 2 

1 1 


1 + e 4 1 + e 2 
1 1 


1 + e 6 1 + e 4 
1 1 

1 + e 8 1 + e 6 
1 1 

1 + e 10 1 + e 8 


1 + e" 


1 


1 + e 

2 

1 


-2 


1 + e' 
1 

1 + e 4 


1 + e° 


1 + e" 


1 


1 + e 

1 


-4 


l + e 

2 

1 


-2 


\ + e A 


l + e A 


l + e 


-8 


l + e 


l + e 


-4 


l + e 

2 

1 


-2 


l + e 2 


l + e 
1 


-10 


l + e 
1 


-8 


l + e' 


l + e 


l + e' 2 
2 


^2 

h 3 

K 

*5 

h 


€ b (L ) — s ^b (* y )| J — 
H^o)- se b {s)\ s=e 
e b (t 0 )-se b (s^ s 


10. Calculate the constant Tj from equation (48). 

11. Obtain the final expression for time-dependent strain by substituting the constants from 
steps 9 and 10 into equation (43). 


D. Determination of Time-Dependent Deflection 

For the one-dimensional beam, the Euler-BernoulH equation for a composite orthotropic 
laminate is expressed as: 


d 2 w 

dx 2 


D*\(t)M x 


(58) 
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For a cantilever beam of constant cross section and thickness, the bending stiffness parameter 

Z>i t (r) is independent of x, and integrating equation (58) twice yields the expression for the deflec- 
tion of the beam. In this manner, one obtains: 


D'u U)P. L 3 

w = . (59) 

3 

The procedure to calculate the time-dependent deflection now follows the one described in 
section IV.C with the following difference. 

After step 7, the load P, like the moment, is divided by the Laplace parameter. Equation (59) 
is then calculated for the value of the Laplace parameter. This is repeated for each value of Laplace 
parameters chosen for the analysis. Equation (51) is again solved with the values of deflection used 
instead of the values of strain. The constants a v and h v are again calculated, and a final expression 
for the time-dependent deflection is then obtained. 


E. Determination of Time- Dependent Natural Frequency 


The natural frequency of a one-dimensional composite orthotropic beam can be determined 
from the elastic beam frequency equation by simply making the following substitution into equation 
( 6 ): 


El = 


b 

D*n W ’ 


( 60 ) 


With this substitution, the resulting expression for natural frequency is: 


O) = 3.51601 


b 

D\ i (t)mL 4 


( 61 ) 


Again the steps to obtain the time-dependent natural frequency are the same as in the time-depen- 
dent strain and deflection with one exception. Notice that the expression for natural frequency is 
independent of applied load. Equation (61), once expressed in the Laplace domain, becomes the 
associated elastic expression. In other words, 5J(.s) = eo(.s). The values of the constants a v and h v 
are again calculated, and the final expression for time-dependent natural frequency also takes the 
form of equation (51). 


V. NUMERICAL EXAMPLES 
A. Single-Phase Viscoelastic Beam Problem 

As a first example of the determination of time-dependent strains and deflection of a can- 
tilever beam of a linear viscoelastic material, the case of a single-phase material will be investi- 
gated. This problem is by no means a new one. It is presented here to demonstrate a simple appli- 
cation of the “correspondence principle” and its usefulness to the design engineer. 
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Figures 2 and 3 show plots of the relaxation modulus and Poisson’s ratio, respectively, for 
the material described in equations (11) and (14). The constants for both equations are defined in 
section III.A. The material presented here is a fictitious one, but the curves represent typical ones 
for actual viscoelastic materials. 

Appendix A. contains the computer program “RELAX” which identifies the step-by-step 
procedure to determine the time-dependent strain and deflection of a linear viscoelastic cantilever 
beam of constant cross section. For this example, the applied load is 25 lb at the tip (free end) of the 
beam. The length of the beam is 29.25 inches, the width is 5 inches, and the thickness is 1 inch. The 
program produces the following results for equations (26) and (29): 

Aj =-0.00514282 
A 2 = -986.113154 
A 3 =-9.859137 

M 0 = 5.555555E-6 
Mj = - 2.698453E-6 
M 2 =- 3.911196E-8 
M 3 = - 4. 023282 E- 8 • 

The plots of strain and deflection from this run of “RELAX” are shown in figures 4 and 5, respec- 
tively. “RELAX” also calculates the deflected shape of the beam as a function of time. This shows 
how the beam relaxes with time, thus yielding a deflection that increases as time passes. The points 
in time selected for the plot are t = 0 h, t = 10 h, t = 100 h, and t = 1,000 h. This is plotted in figure 6. 


B. Two-Phase Viscoelastic Beam Problem 

Many composite materials have, as constituents, a time-dependent component (polymer 
matrix) and fibers or particulates that are made from materials that have properties that are rela- 
tively insensitive to creep or relaxation (graphite, boron, silicone-carbide, etc.). In these cases, the 
problem of determining the time-dependent response to applied loads becomes more involved, and in 
many cases, closed-form solutions are not available. The intent of this section is to use the corre- 
spondence principle and apply the collocation method 13 to a laminated composite cantilever beam 
loaded at the free end to determine the strain, deflection, and natural frequency as a function of time. 
The application of the procedure to different problems has been well documented. 8 9 The approach is 
given in section IV and is applied and explained in a step-by-step manner in the computer program 
“VISCOBM” found in appendix B. 
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Figure 5. Time-dependent deflection for a single-phase viscoelastic beam 






3.00 



x-coordinate of beam (in) 


Figure 6. Deflected shape of a single-phase viscoelastic cantilever beam at various times. 

One will notice, when using the collocation method, that the accurate determination of the 
initial response ( t = 0) is critical to the proper solution of the viscoelastic problem. For this purpose, 
a computer program “VISTO” is included in appendix C. This program is nothing more than 
“VISCOBM” at time t = 0. It is important to point out that it is not necessary to create a separate 
program for this condition since running “VISCOBM” at t = 0 would accomplish the same results. It 
is done this way here as a means of simply demonstrating the. method. 

Figure 7 shows the calculated time-dependent strain for various ply lay-up angles of a four- 
ply symmetric beam. It is interesting to see how the lay-up angle plays a significant role in the long- 
term strain values. Figure 8 shows the same trends for deflection. Figure 9 illustrates that increas- 
ing the fiber content of the composite not only increases the stiffness of the beam but it also reduces 
the long-term effects of relaxation in the beam. Figure 10 shows the significant drop in natural 
frequency of the beam as a function of time. 


VI. LIMITATIONS OF THE CORRESPONDENCE PRINCIPLE 


A very important limitation to the application of the correspondence principle must be 
explained. The knowledge of the state of the boundary conditions must be known in order to apply 
this method effectively. 15 

If the interface between the surfaces where the stress is prescribed and where the displace- 
ments are prescribed changes with time, the correspondence principle is not applicable. The condi- 
tions of each surface, however, can be time-dependent. The illustration in figure 1 1 shows how the 
shaded area (interface boundary area) in a viscoelastic medium changes as a function of time for a 
spherical indentor. For a cylindrical indentor, the interface boundary does not change. An example of 
a changing interface boundary is the increasing inner diameter of a solid propellant motor as it is 
consumed during operation. 
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Figure 7. Time-dependent bending strain of outer ply in a two-phase viscoelastic cantilever beam 

for various laminate configurations. 
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Figure 8. Normalized time-dependent deflection of a two-phase viscoelastic cantilever beam 

for various laminate configurations. 
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vn. CONCLUSIONS 


The elastic- viscoelastic correspondence principle is a very effective and straightforward tool 
for determining the time-dependent response of viscoelastic structures to applied loads. The design 
engineer can use this knowledge to optimize parameters such as ply lay-up angle, fiber volume 
content, and configuration variables. This can help to minimize or maximize the effects of creep or 
relaxation depending on the functional use of the hardware designed. The examples presented offer 
the design engineer the capability to understand the effect that a viscoelastic material can have on 
the performance of structures. 

Although not included here, temperature and exposure to environmental effects (ultraviolet 
light, ozone, atomic oxygen, etc.) can play an important role in the degradation of the material 
properties. These effects should be included in the development of die basic material properties such 
as the relaxation modulus, creep compliance, and Poisson’s ratio. 

As composite materials become more widely used in aerospace applications, the effects of 
long-term exposure to the space environment can become a significant factor in the geometry and 
specifications of the hardware. With the method presented here, the design engineer has the tools 
necessary to alter initial designs which have not taken into consideration the time-dependency 
factor. This will provide for more structurally sound and efficient structures. 
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APPENDIX A 


COMPUTER PROGRAM “RELAX” 


on non 


C*********************************************************************** 

C RELAX 

C***************************************** ******************** ********** 

C PROGRAM TO DETERMINE THE TIME DEPENDENT STRAIN AND 

C DEFLECTIONS OF A SINGLE PHASE VISCOELASTIC 

C CANTILEVER BEAM 

C*************** ********************************************** ********** 

IMPLICIT DOUBLE PRECISION (A-H, O-Z) 

DIMENSION XL (4) ,XLMD(3) , RI ( 3 ) , FA (3 ) , GA (3 ) , XM(3 ) 

OPEN (UNIT=7 , FILE= * C : \FORTRAN\ RELAX . DAT 1 ) 

OPEN(UNIT=8, FILE= 1 C : \FORTRAN\STRAIN.DAT ' ) 

C 

C IDENTIFY BEAM GEOMETRY AND APPLIED LOAD 

C 

P = 25.00 
XLTH = 29.25 
WIDE =5.00 
THK =1.00 

c 

C CALCULATE MOMENT OF INERTIA AND SPECIFIC CONSTANTS 

C 

XINRT = {WIDE* (THK**3) ) /12 . 

F1XZ = 12 *P*XLTH* (THK/2 ) / {WIDE* (THK**3) ) 

Q = 4 . *P* (XLTH**3 ) / (WIDE* {THK**3 ) ) 

T = 0.0001 


IDENTIFY MATERIAL CONSTANTS FOR RELAXATION MODULUS FUNCTION 


A = 180000. 
B1 = 5000. 

B2 = 5000. 

B3 = 170000. 
GAMA1 = 1000 
GAMA2 = 10. 
GAMA3 = .01 


C 


C 

C 

C 

C 


CALCULATE RELAXATION MODULUS 


DO 20 I = 1,1301 

E = A + Bl*EXP{-GAMAl*T) + B2*EXP{-GAMA2*T) + B3 * EXP ( -GAMA3 *T ) 
WRITE ( 7 , * ) T, E 
T = T + 1 
20 CONTINUE 


CALCULATE CONSTANTS FOR CALCULATION OF STRAIN AND 
DEFLECTION IN THE LAPLACE DOMAIN 


FI = GAMA1 + GAMA2 + GAMA3 

F2 = GAMA1 *GAMA2 + GAMA2 *GAMA3 + GAMA3 *GAMA1 

F3 = GAMA1 *GAMA2 *GAMA3 

F4 = GAMA2 + GAMA3 

F5 = GAMA2 *GAMA3 

F6 = GAMA1 + GAMA3 

F7 = GAMA1 *GAMA3 
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F8 = GAMA1 + GAMA2 

F9 = GAMA1 *GAMA2 

XL ( 1 ) = A + B1 + B2 + B3 

XL ( 2 ) = A*F1 + B1 *F4 + B2*F6 + B3*F8 

XL (3 ) = A*F2 + B1*F5 + B2*F7 + B3*F9 

XL {4} = A*F3 

C 

C OBTAIN THE ROOTS OF THE Ga(s) TERM IN THE DENOMINATOR OF EQUATION (20) 

C 

CALL CUBIC (XL, XLMD, RI ) 

WRITE (7, lO)XLMD(l) ,XLMD(2) , XLMD (3) 

10 FORMAT (' THE REAL ROOTS OF Ga(s) ARE\3E12.5) 

WRITE (7, ll)RI(l) ,RI(2) ,RI(3) 

11 FORMAT {' THE IMAGINARY ROOTS OF Ga(s) ARE',3E12.5) 

C 

CALCULATE THE "M“ CONSTANTS OF EQUATION (24) BY HEAVISIDE'S 
PARTIAL FRACTION EXPANSION THEOREM 


XMO = F3 /XL ( 4 ) 

XXMO = F3 /XXL ( 4 ) 

WRITE (7, 12) XMO 

12 FORMAT (' THE CONSTANT Mo IS\E12.5) 

DO 21 I = 1,3 

FA ( I ) = XLMD ( I ) **3 + Fl *XLMD ( I ) **2+F2 *XLMD ( I ) +F3 

GA ( I ) = 4 . *XL ( 1 ) *XLMD { I ) * *3 + 3 . *XL ( 2 ) *XLMD ( I ) * *2+2 . *XL { 3 ) *XLMD (I) + 
*XL ( 4 ) 

XM ( I ) = FA ( I ) /GA ( I ) 

WRITE (7 , 13) I , XM ( I ) 

13 FORMAT {' M ' , I 1 , 1 =',E12.5) 

21 CONTINUE 


CALCULATE THE TIME DEPENDENT STRAIN PER EQUATION (26) 


T = .0001 

DO 22 I = 1,2500,5 

STRN = F1XZ* (XM ( 1 ) *EXP (XLMD ( 1 ) *T) +XM ( 2 ) *EXP (XLMD{2 ) *T) + 
*XM ( 3 ) *EXP (XLMD (3 ) *T) +XM0) 


CALCULATE THE TIME DEPENDENT DEFLECTION PER EQUATION (29) 


DEFL = Q* (XM ( 1 ) *EXP (XLMD ( 1 ) *T) +XM ( 2 ) *EXP (XLMD (2 ) *T) + 
*XM ( 3 ) *EXP (XLMD ( 3 ) *T) +XM0 ) 

WRITE (8, 30) T, STRN, DEFL 
30 FORMAT (11E12 .5) 

T = T + 5 
22 CONTINUE 


CALCULATE THE DEFLECTED SHAPE OF THE CANTILEVER BEAM 
AT SEVERAL DIFFERENT TIMES (0 hrs, 10 hrs, 100 hrs & 1000 hrs) 


X = 0.0 

DO 23 I = 1,21 

Q1 = ( P*XLTH* (X* *2 ) /2 . - P* (X**3 ) /6 . ) /XINRT 
Q1X = ( P*XLTH* (X**2 ) /2 . - P*(X**3)/6.) 


C 

C T = 0.0 hrs. 

C 

T=0 

DO = Q1 * (XM { 1 ) *EXP (XLMD ( 1 ) *T ) +XM ( 2 ) *EXP (XLMD ( 2 ) *T ) + 

*XM (3 ) *EXP (XLMD (3 ) *T) +XM0 } 

c 

C T = 10.0 hrs. 

C 

T=10 

DIO = Ql* (XM ( 1) *EXP(XLMD( 1) *T) +XM(2 ) *EXP (XLMD (2 ) *T) + 
*XM (3 ) * EXP (XLMD ( 3 ) *T) +XM0 } 

C 

C T = 100.0 hrs. 

C 

T=100 

D100 = Ql* (XM { 1 ) *EXP (XLMD { 1 ) *T) +XM ( 2 ) *EXP {XLMD (2 ) *T) + 
*XM(3 ) *EXP ( XLMD ( 3 ) *T ) +XM0 ) 

c 

C T = 1000.0 hrs. 

C 

T=1000 

D1000 = Ql * { XM ( 1 ) *EXP (XLMD ( 1 ) *T ) +XM ( 2 ) *EXP (XLMD { 2 ) *T) + 
*XM(3 ) *EXP (XLMD (3 ) *T) +XM0 ) 

WRITE ( 9, 30) X, DO f D10 # D10 0 , D10 00 
X = X + 1.4625 
23 CONTINUE 
STOP 
END 


28 


APPENDIX B 


COMPUTER PROGRAM “VISCOBM” 


n o 


*************************** ******************************************** 

PROGRAM TO DETERMINE THE TIME DEPENDENT STRAINS , DEFLECTION 
C AND NATURAL FREQUENCY OF A COMPOSITE BEAM SUBJECTED TO BENDING LOADS 

C* ********************************************************************** 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION ALPHA (4) 

DIMENSION E ( 4 , 22 ) 

DIMENSION XNU (4,21) , G { 4 , 12) 

DIMENSION ZT { 5 ) 

DIMENSION D(6, 6) / EPSX{4) ,EPSY{4) , EPSXY{4) 

DIMENSION SIGX ( 4 ) ,THK(4) ,SIGY(4) / SIGXY(4) 

DIMENSION SL ( 6 ) , EMS { 6 ) , PMS ( 6 ) 

DIMENSION STS1 (200) , STS2 (200) , STS3 (200) ,STS4 (200) 

DIMENSION HI (6) , H2 (6) , H3 (6) , H4 (6) 

DIMENSION F1<6),F2(6) ,F3(6) ,F4{6) ,ALAP(6,6) 

DIMENSION ALAP2 (6, 6) , ALAP3 (6,6) , ALAP4 (6,6) 

DIMENSION QR1 1(4) , QR22 (4) ,QR12 (4) , QR66(4) , QR26 (4) , QR16 (4 ) 
DIMENSION SRI 1(4) , SR22 (4) , SR12 (4) , SR66 (4) ,SR26(4) ,SR16(4) 
DIMENSION Sll (4) ,S22 (4) # S12 (4) ,S66(4) 

DIMENSION Qll (4) ,Q22(4) ,Q12 (4) ,Q66(4) 

DIMENSION XIN (31) , CAL (31) ,XLAMD(4) , F5(6) , ALAP5 ( 6 , 6 ) , H5 ( 6 } 
DIMENSION F6 (6) , ALAP6 (6, 6) , H6 (6) 

C — 

C IDENTIFY OUTPUT FILES 

C 

OPEN (UNIT=12 , FILE= ' C : \FORTRAN\ ACTUAL . DAT ' ) 

OPEN (UNIT=13 , FILE= 1 C : \FORTRAN\ NORMAL . DAT ' ) 

c — 

C INPUT WIDTH AND LENGTH OF BEAM 


WIDTH 

= 

4. 

XLB = 

29 

.25 

INPUT 

VALUES OF THE LAPLACE PARAMETER AT DISCRETE POINTS 

SL(1) 

— 

.0067379 

SL ( 2 ) 

- 

.0497871 

SL(3) 

= 

.3678794 

SL ( 4 ) 

= 

2.7182818 

SL (5) 

= 

20.0855369 

SL(6) 

= 

148.4131591 


C 

q ************************************************* 

C * IDENTIFY LAMINA MATERIAL PROPERTIES * 

Q ************************************************* 

c 

C INPUT VOLUME FRACTION FOR MATRIX AND FIBERS (VM & VF) 

C INPUT BULK MODULUS OF ELASTICITY FOR MATRIX <XK) 

C INPUT POISSON'S RATIO AND YOUNG 1 S MODULUS FOR FIBERS (XNUF & EF) 

C INPUT TOTAL NUMBER OF PLIES ( NLAY ) 

C INPUT DENSITY OF A PLY (RHO) 

C 

VM = .38 
VF = .62 
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XK = 400000. 
EF = 33.33E6 
XNUF = .2 
NLAY = 4 
PI = 3.14159 
RHO = .057 


CALCULATE SHEAR MODULUS FOR FIBERS 


GF = EF/ (2* (1+XNUF) ) 


INPUT PREVIOUSLY CALCULATED VALUE OF PLY BENDING STRAINS, 
DEFLECTION AND NATURAL FREQUENCY AT INITIAL TIME (T=0.0) 

(FROM VI ST 0. FOR) 


THIS CASE IS 4 PLY [3 0/-30] symmetric 


EPS01 

= 

. 00140303 

EPS02 

= 

.000701517 

EPS 03 

= 

. 00000000 

EPS 04 

= 

-.000701517 

DEFT0 

= 

3.201 

FREQ0 

= 

10 . 1 


INPUT PRONY COEFFICIENTS FOR RELAXATION MODULUS 
C AND TIME DEPENDENT POISSON'S RATIO 

C 

C RELAXATION MODULUS COEFFICIENTS 

C 

A1 = 180000. 

B1 = 5000. 

B2 = 5000. 

B3 = 170000. 

C 

C POISSON'S RATIO COEFFICIENTS 

c 

A13 = .425 
B13 = -.0020833 
B23 = -.0020833 
B33 = -.070833 

c 

C CALCULATE THE RELAXATION MODULUS AND TIME DEPENDENT 

C POISSON'S RATIO FOR THE COMPOSITE SYSTEM MATRIX MATERIAL 

C IN THE LAPLACE DOMAIN 

C 

C EMS(I) = RELAXATION MODULUS FOR MATRIX 

C PMS(I) = POISSON'S RATIO FOR MATRIX 

C SL ( I ) = LAPLACE PARAMETER 

c 

C 

DO 1000 LAP =1,6 

EMS ( LAP ) = A1 + SL (LAP) *B1/ (SL (LAP) +1000 . ) + SL(LAP) *B2/ (SL (LAP) 
1 10.) + SL(LAP)*B3/(SL(LAP)+.01) 

PMS (LAP) = A13 + SL(LAP) *B13/ (SL(LAP) +1000 . ) + 


o n 


1 SL (LAP) *B23/ (SL (LAP) +10 . ) + SL (LAP) *B33/ (SL (LAP) + . 01 ) 


CALCULATE SHEAR MODULUS FOR THE MATRIX 

C 

GM = EMS (LAP) / (2* ( 1+PMS (LAP) ) ) 

C 

C CALCULATE LONGITUDINAL AND TRANSVERSE PROPERTIES OF A 

C UNIDIRECTIONAL LAMINA (USING MICROMECHANICS APPROACH) 

C 

DO 10 I = 1 , NLAY 

E(I, 11) = EF*VF + EMS(LAP)*VM 

E(I, 22) = EF*EMS (LAP) / (VM*EF + VF*EMS(LAP)) | 

GNUM = GF* ( PI + 4 *VF ) + GM*(PI-4*VF) 

GDEN = GF* ( PI-4*VF) + GM*(PI+4*VF) 

GN = GNUM/GDEN 
G121 = ( (4-PI)+PI*GN) /4 
G122 = 4*GN/ ( (4-PI) *GN+PI) 

G ( I , 12 ) = (GM/2 ) * (G12 1 + G122 ) 

XNU ( I , 12) = XNUF*VF + PMS(LAP)*VM 
XNU (1,21) = XNU (I, 12) *E(I,22)/E(I, 11) 

10 CONTINUE 

c 

C INPUT THE ANGULAR ORIENTATION OF EACH LAMINA WITH RESPECT TO 

C THE LONGITUDINAL AXIS OF THE LAMINATE 

C 

ALPHA ( 1 ) = 30.0 

ALPHA ( 2 ) = -30.0 

ALPHA ( 3 ) = -30.0 

ALPHA ( 4 ) = 30.0 

C 

C CALCULATE THE REDUCED STIFFNESS COEFFICIENTS FOR EACH LAMINA 

C 

DO 11 I = 1 , NLAY 

Qll(I) = E ( I , 1 1 ) / ( 1 . -XNU ( I , 12 ) *XNU (1,21) ) 

Q12 { I ) = XNU (I,12)*E(I,22)/(1. -XNU ( I , 12 ) *XNU (1,21) ) 

Q22 ( I ) = E(I, 22) / (1 .-XNU(I, 12) *XNU ( 1 , 21) ) 

Q66 ( I ) = G(I, 12) 

C 

C CALCULATE THE COMPLIANCE MATRIX COEFFICIENTS FOR EACH LAMINA 

C 

Sll(I) = 1/E { I , 11 ) 

S12 { I ) = -XNU ( I, 12 ) /E( I, 11) 

S22 ( I ) = 1/E { I , 22 ) 

S66 (I) = 1/G { I , 12 ) 

c 

C CALCULATE THE TRANSFORMED REDUCED STIFFNESSES FOR EACH LAMINA 



ALPHA ( I ) = ALPHA (I) *PI/180 . 

QRll(I) = Qll ( I ) * (COS (ALPHA ( I ) ) ) * *4 

1 +2 . * (Q12 ( I ) +2 . *Q66 ( I ) ) * ( (SIN (ALPHA ( I ) ) )**2) * 

2 { (COS { ALPHA ( I ) ) ) **2)+Q22 (I) * (SIN(ALPHA(I) ) ) **4 
C 

QR12 ( I ) = (Qll ( I ) +Q22 ( I ) -4 . *Q66 ( I ) ) * ( (SIN ( ALPHA ( I) ) ) **2 ) * 

1 ( (COS (ALPHA(I) ) ) **2 ) +Q12 ( I ) * ( (SIN (ALPHA ( I ) ) ) **4+ 
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2 ( COS ( ALPHA ( I ) ) ) * * 4 ) 

C 

QR 16(1) = (Qll ( I ) -Q12 ( I) -2 . *Q66 (I) ) * ( (COS (ALPHA ( I ) ) ) **3 ) * 

1 SIN ( ALPHA ( I ) ) + (Q12 { I ) -Q2 2 ( I ) +2 . *Q66 ( I ) ) *COS (ALPHA ( I) ) * 

2 (SIN (ALPHA (I) ) ) **3 
C 

QR22 ( I ) = Qll (I) * (SIN (ALPHA (I) ) ) **4 

1 +2 . *(Q12 (I) +2 . *Q66(I) ) * ( (SIN (ALPHA (I) ) ) **2) * 

2 ( (COS (ALPHA ( I) ) ) **2)+Q22 (I) * (COS (ALPHA ( I ) ) ) **4 
C 

QR2 6(1) = (Qll (I) -Q12 (I) -2 . *Q66(I) ) * ( (SIN (ALPHA ( I ) ) ) **3) * 

1 COS (ALPHA (I) ) + (Q12 (I) -Q22 (I) +2 . *Q66 (I) ) *SIN ( ALPHA ( I ) ) * 

2 ( COS ( ALPHA ( I ) ) ) * * 3 
C 

QR66(I) = (Qll (I) +Q22 (I) -2 . *Q12 (I) -2 . *Q66 ( I) ) * ( (SIN (ALPHA ( I ) ) )**2) 

1 *( (COS (ALPHA ( I ) ) )**2)+Q66(I)*( (SIN (ALPHA { I ) ) )**4+ 

2 ( COS ( ALPHA ( I ) ) ) * * 4 ) 

C 

C CALCULATE THE TRANSFORMED COMPLIANCE COEFFICIENTS FOR EACH LAMINA 

c 

SR 11 ( I ) = Sll ( I) * (COS ( ALPHA { I ) ) ) * * 4 

1 + (2 . *S12 ( I ) + S66 ( I ) ) * ( (SIN (ALPHA ( I) ) ) **2 ) * 

2 ( (COS (ALPHA (I) ) ) **2)+S22 (I) * (SIN (ALPHA { I ) ) ) **4 
C 

SR12 ( I ) = (Sll (I)+S22 (I) -S66 (I) ) * { (SIN (ALPHA { I ) ) ) **2 ) * 

1 ( (COS (ALPHA (I) ) ) **2 ) +S12 ( I ) * ( (SIN (ALPHA ( I ) ) )**4+ 

2 (COS (ALPHA ( I ) } ) * *4 ) 

C 

SRI 6 ( I ) = (2 . *S11 (I) -2 . *S12 (I) -S66 (I) ) *( (COS (ALPHA ( I ) ) ) **3) * 

1 SIN (ALPHA (I) ) - (2*S22 (I) -2*S12 (I) -S66(I) ) *COS ( ALPHA ( I ) 

2 ) *(SIN(ALPHA(I) ) ) **3 
C 

SR22 ( I ) = S11(I)*(SIN (ALPHA ( I ) ) ) * *4 

1 + (2*S12 (I)+S66 (I) ) * ( (SIN (ALPHA ( I ) ) ) **2) * 

2 ( ( COS ( ALPHA ( I ) ) ) **2)+S22 (I) * (COS (ALPHA ( I ) ) ) **4 
C 

SR26 { I ) = (2*Sll (I) -2*S12 (I) -S66 (I) ) * ( (SIN (ALPHA { I ) ) ) **3 ) * 

1 COS ( ALPHA { I ) ) - (2*S22 (I) -2*S12 (I) -S66 ( I ) ) *SIN ( ALPHA ( I ) 

2 ) * ( COS ( ALPHA ( I ) ) ) * * 3 
C 

SR66 ( I ) = 2* (2*S11 (I)+2*S22 (I) -4 . *S12 (I) -S66 (I) ) * { (SIN ( ALPHA ( I ) ) 

1 ) **2) * ( (COS (ALPHA (I) ) } **2)+S66(I) *( (SIN (ALPHA ( I ) ) ) **4+ 

2 ( COS ( ALPHA ( I ) ) ) * * 4 ) 

11 CONTINUE 

Q ******************************************************* 

C * RELATE EACH LAMINA THICKNESS TO ITS LOCATION WITHIN * 

C * THE LAMINATE (Z COORDINATE IS 0 AT THE MIDSURFACE) * 

C 

C INPUT LAMINA THICKNESSES 

C 

DO 998 J = 1 , NLAY 
THK(J) = .0625 
998 CONTINUE 
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c 

TOTTH =0.0 
DO 25 J = 1 , NLAY 
TOTTH = TOTTH + THK(J) 

25 CONTINUE 

C 

C CALCULATE THE TOTAL WEIGHT OF THE BEAM 

C 

TOTWT = WIDTH *TOTTH*XLB*RHO 

C 

HAFTH = TOTTH/ 2 . 

ZT (NLAY) = HAFTH 
DO 24 J = 1 , NLAY 
LAM = NLAY - J + 1 
ZT(NLAY-J) = ZT (LAM) - THK(LAM) 

24 CONTINUE 

DO 40 J = 1 , NLAY + 1 
40 CONTINUE 

C 

C CALCULATE BENDING STIFFNESSES FOR THE ENTIRE LAMINATE 

C 

Dll = 0.0 

D12 =0.0 

D16 =0.0 

D22 =0.0 

D26 =0.0 

D66 =0.0 

DO 35 J = 1 , NLAY 

L = J + 1 

Dll = Dll + QR11 ( J) * ( (ZT (L-l ) **3 . ) - (ZT (L-2 ) **3 . ) ) /3 . 

D12 = D12 + QR12 ( J) * ( (ZT(L-1)**3.)- (ZT (L-2 ) **3 . ) ) /3 . 

D16 = D16 + QR16 (J)*((ZT(L-l)**3.)-(ZT(L-2)**3.))/3. 

D22 = D22 + QR22 (J) * ( (ZT(L-l) **3 . ) - (ZT (L-2 ) **3 . ) ) /3 . 

D26 = D26 + QR26 (J)*-((ZT(L-l)**3.)-(ZT(L-2)**3.))/3. 

D66 = D66 + QR66 ( J) * ( (ZT (L-l ) * *3 . ) - (ZT (L-2 ) **3 . ) ) /3 . 

35 CONTINUE 

C — 


C CALCULATE THE DETERMINANT OF THE STIFFNESS MATRIX j 

c— — — — - — ----- — — — } 

DETER = Dll*D22*D66 + D12*D26*D16 + D16*D12*D26 j 

1 - D16*D16*D22 - D12*D12*D66 - D26*D26*D11 f 

C — ■ "i 


C CALCULATE COEFFICIENTS FOR THE COFACTORS OF THE STIFFNESS MATRIX 

C 

DllC = D22*D66 - D26*D26 
D12C = - (D12*D66 - D16*D26) 

D16C = D12*D26 - D16*D22 

C — 

C CALCULATE THE COEFFICIENTS OF THE INVERSE STIFFNESS MATRIX 

C 

D11S = D11C/DETER 
D12S = D12C/DETER 
D16S = D16C/ DETER 
WRITE ( 13 , *) ’ D11S = ' , D11S 
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IDENTIFY APPPLIED MOMENTS 


THE LOAD MUST BE EXPRESSED IN THE LAPLACE DOMAIN 
LOAD ( S ) = LOAD/SL (LAP) [for a CONSTANT load] 


XMAP = -292.5/WIDTH 
XMAP = XMAP/SL (LAP) 
PLOAD = 10 . /SL (LAP) 


*************** PERFORM DEFLECTION CALCULATIONS ****************** 


CALCULATE THE CURVATURE AT ANY POINT ALONG THE LENGTH OF THE BEAM 


INM1 =30 
XIN{1) =0.0 
DO 805 K = 2,31 
XIN(K) = XIN (K-l ) +XLB/ INM1 
805 CONTINUE 

DO 804 NI = 1, INM1+1 

CAL(NI) = D11S*PL0AD* (XLB-XIN (NI) ) /WIDTH 
804 CONTINUE 


PERFORM A POLYNOMIAL CURVE FIT FOR THE CURVATURE AT 
ANY POINT ALONG THE LENGTH 


CALL POLYCF ( X IN , C AL , 3 1 , 3 , XL AMD , NERR ) 


CALCULATE THE TIP DEFLECTION OF THE BEAM BY INTEGRATING THE 
CURVATURE TWICE WITH RESPECT TO THE LENGTH COORDINATE 


DEF =0.0 
DO 111 J = 1,4 

DEF = DEF + XLAMD(J) * (XLB** (J+l) ) / (J* (J+l) ) 
111 CONTINUE 


CALCULATE THE ASSOCIATED ELASTIC NATURAL FREQUENCY OF THE BEAM 


FREQ = 3 . 51 601 *SQRT ( (WIDTH*386 . 4 ) / (D1 1S*T0TWT* (XLB* *4 ) ) ) 


IDENTIFY THE DEFLECTION AS THE ASSOCIATED ELASTIC 
SOLUTION MULTIPLIED BY THE LAPLACE PARAMETER 


DEF = DEF*SL ( LAP) 


CALCULATE THE STRESSES FOR EACH LAMINA 


DO 36 J = 1 , NLAY 
L = J + 1 

SIGX(J) = ZT (L-2 ) *XMAP* (QR11 ( J) *D11S +QR12 { J) *D12S +QR16 ( J) *D16S) 
SIGY(J) = ZT (L-2 ) *XMAP* (QR12 ( J) *D11S +QR22 ( J) *D12S +QR26 ( J) *D16S) 
SIGXY(J) = ZT (L-2 ) *XMAP* (QR16 ( J) *D1 IS +QR26 ( J) *D12S +QR66 ( J) *D16S) 
36 CONTINUE 
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c 

C CALCULATE THE TIME DEPENDENT STRAINS FOR EACH LAMINA 

C 

DO 37 J = 1,NLAY 

EPSX(J) = SR11 { J) *SIGX (J) + SR12 { J) *SIGY { J) +SR16 (J) *SIGXY ( J) 

EPSY(J) = SR 12 { J) *SIGX ( J) + SR2 2 ( J) *SIGY ( J) +SR26 ( J) *SIGXY ( J) 
EPSXY(J) = SRI 6 ( J) *SIGX { J) +SR26 { J) *SIGY { J) +SR66 (J) *SIGXY(J) 

37 CONTINUE 

C 

C IDENTIFY STRAINS FOR EACH LAMINA AS THE ASSOCIATED ELASTIC 

C STRAINS MULTIPLIED BY THE LAPLACE PARAMETER 

C 

DO 38 J = 1,NLAY 
EPSX(J) = EPSX(J) *SL(LAP) 

EPSY(J) = EPSY(J) *SL(LAP) 

EPSXY(J) = EPSXY ( J) *SL (LAP) 

38 CONTINUE 

C 

C CALCULATE DIFFERENCE BETWEEN STRAINS, DEFLECTION AND NATURAL 

C FREQUENCY AT TIME ZERO AND AT THE DISCRETE VALUES 

C OF THE LAPLACE PARAMETER FOR INVERSION BACK INTO 

C THE TIME DOMAIN 

C — “r-r 

C NOTE THAT ONLY THE STRAINS IN THE X DIRECTION ARE USED 

C FOR THIS PROBLEM 

C 

FI (LAP) = EPS01 - EPSX(l) 

F2 (LAP) = EPS02 - EPSX(2) 

F3 (LAP) = EPS03 - EPSX(3) 

F4 ( LAP) = EPS04 - EPSX(4) 

F5 ( LAP) = DEFTO - DEF 

F6 ( LAP) = FREQO - FREQ 

WRITE ( 13 , 9 02 ) LAP, FREQO , FREQ, F6 (LAP) 

902 FORMAT ( *LAP= ',12, 1 FREQ0= ' , E12 . 5 , ' FREQ= ' , E12 . 5, 1 F6 (LAP) = 1 , E12 . 5 ) 
1000 CONTINUE 


PERFORM LAPLACE TRANSFORM INVERSION AND SOLVE FOR THE CONSTANTS 
C NEEDED FOR THE CALCULATIONS OF TIME DEPENDENT STRAINS 

C 

DO 1001 LAPW =1,6 
DO 1001 LAPV =1,6 
ALPWV = EXP (2* (LAPW-LAPV) ) 

ALAP( LAPW, LAPV) = 1/(1 + ALPWV) 

ALAP2 (LAPW, LAPV) = 1/(1 + ALPWV) 

ALAP3 (LAPW, LAPV) = 1/(1 + ALPWV) 

ALAP4 (LAPW, LAPV) = 1/(1 + ALPWV) 

ALAP5 (LAPW, LAPV) = 1/(1 + ALPWV) 

ALAP6 (LAPW, LAPV) = 1/(1 + ALPWV) 

1001 CONTINUE 

C 

C SOLVE SIMULTANEOUS EQUATIONS 

C — — — 

CALL MATINV ( ALAP, 6 , FI , HI ) 

CALL MATINV (ALAP2, 6, F2,H2) 
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CALL MATINV ( ALAP3 , 6 , F3 , H3 ) 
CALL MATINV {ALAP4 , 6 , F4 , H4 ) 
CALL MATINV (ALAP5, 6, F5, H5) 
CALL MATINV ( ALAP6 , 6 , F6 , H6 ) 


: EVALUATE GAMMA 

CONSTANTS FROM INITIAL STRAIN CONDITIONS 

SUMSTl 

o 

o 

II 


SUMST2 

= 0.0 


SUMST3 

o 

o 

II 


SUMST4 

= 0.0 


SUMST5 

o 

o 

II 


SUMST6 

II 

o 

o 


DO 1002 

LAP =1,6 

SUMSTl 

= SUMSTl 

+ HI (LAP) 

SUMST2 

= SUMST2 

+ H2 (LAP) 

SUMST3 

= SUMST3 

+ H3 (LAP) 

SUMST4 

= SUMST4 

+ H4 (LAP) 

SUMST5 

= SUMST5 

+ H5 (LAP) 

SUMST6 

= SUMST6 

+ H6 (LAP) 

1002 CONTINUE 


GAMA1 = 

EPS01 - 

SUMSTl 

GAMA2 = 

EPS02 - 

SUMST2 

GAMA3 = 

EPS03 - 

SUMST3 

GAMA 4 = 

EPS04 - 

SUMST4 

GAMAS = 

DEFT0 - 

SUMST5 

GAMA 6 = 

FREQ0 - 

SUMST6 


CALCULATE TIME DEPENDENT STRAINS , DEFLECTION 
AND NATURAL FREQUENCY 

LTIME = 

200 


TIME = 

.0001 



DO 1003 I = 1 , LTIME 
STS1 ( I ) = 0.0 
STS2 ( I ) =0.0 
STS3 ( I ) =0.0 
STS 4 { I ) =0.0 
DEFLEC =0.0 
FREQCY = 0.0 
DO 1004 LAP =1,6 

STSl(I) = STSl(I) + HI (LAP) *EXP( -SL(LAP) *TIME) 
STS2 { I ) = STS2 { I ) + H2 (LAP) *EXP ( -SL (LAP) *TIME) 
STS3 ( I ) = STS3 ( I ) + H3 ( LAP) *EXP ( -SL (LAP) *TIME) 
STS4 ( I ) = STS4 { I ) + H4 (LAP) *EXP ( -SL (LAP) *TIME) 
DEFLEC = DEFLEC + H5 (LAP) *EXP { -SL (LAP) *TIME) 
FREQCY = FREQCY + H6 ( LAP) *EXP ( -SL ( LAP) *TIME) 
1004 CONTINUE 

STS1(I) = GAMA1 + STSl(I) 

STS2 ( I ) = GAMA2 + STS2(I) 

STS3 ( I ) = GAMA3 + STS3(I) 

STS4 ( I ) = GAMA4 + STS4 { I ) 

DEFLEC = GAMA5 + DEFLEC 
FREQCY = GAMA 6 + FREQCY 
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WRITE ( 12 , 1006) I , TIME, STS1 ( I ) , STS2 { I ) , STS3 (I) , STS4 ( I ) , DEFLEC, FREQCY 
TIME = TIME* 1 . 1 
1003 CONTINUE 

C 

C NORMALIZE THE STRAINS WITH RESPECT TO THE MAXIMUM 

C FOR PURPOSES OF GRAPHING THE RESULTS 

C (TMAG = MAGNIFICATION FACTOR) 

C — 

TIME = ,0001 
TMAG = 100, 

DO 285 I = 1 # LTIME 

STS1 (I) =TMAG*STS1 (I) /STS1 (LTIME) 

STS2 (I) =TMAG*STS2 (I) /STS2 (LTIME) 

C** ***** ip the STFAIN is zero (mid surface) DO NOT NORMALIZE*********** 
STS3 ( I ) IS THE STRAIN AT THE MIDDLE SURFACE 
********************************************************************** 

STS4 (I) =TMAG*STS4 (I) /STS4 (LTIME) 

WRITE (13, 1006) I , TIME, STS1 ( I ) , STS2 (I) , STS3 (I) ,STS4(I) 

1006 FORMAT ( 15 , 7E12 .5) 

TIME = TIME*1 . 1 
285 CONTINUE 
STOP 
END 
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APPENDIX C 


COMPUTER PROGRAM “VISTO” 



n a 


c*********************************************************************** 

C* PROGRAM TO DETERMINE THE INITIAL (T=0.0) STRAINS AND DEFLECTION * 
C* OF A VISCOELASTIC COMPOSITE BEAM SUBJECTED TO BENDING LOADS * 

Q* ********************************************************************** 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION ALPHA (24) 

DIMENSION E ( 24 , 22 ) 

DIMENSION XNU(24, 21) ,G(24, 12) 

DIMENSION Z (25) , ZT (25) 

DIMENSION D ( 6 , 6) ,EPSX(24) ,EPSY{24) ,EPSXY(24) 

DIMENSION SIGX (24 ) ,THK(24) ,SIGY(24) ,SIGXY(24) 

DIMENSION QR11 (24) , QR2 2 (24) , QR12 (24) ,QR66(24) ,QR26(24) ,QR16(24) 
DIMENSION SRI 1 (24) , SR22 (24) , SR12 (24) , SR66 (24) , SR26 (24) , SRI 6 (24) 
DIMENSION Sll (24) ,S22 (24) ,S12 (24) , S66 ( 24 ) 

DIMENSION Qll (24) ,Q22 (24) ,Ql2 (24) ,Q66(24) 

DIMENSION XIN (31) , CAL (31) , XLAMD ( 4 ) 

OPEN (UNIT=13 , FILE= ' C : \FORTRAN\VISTO . DAT 1 ) 



I FLAG = 0 


C 


C 


C 

C 

C 


C 

C 

C 


IDENTIFY LAMINA MATERIAL PROPERTIES 
AT TIME T = 0 . 0 SECONDS 


PI = 3.14159 
WRITE (6 , *) 1 INPUT 
READ { 5 , * ) VM 
WRITE (6, *) ' INPUT 
READ ( 5 , * ) VF 
WRITE ( 6 , * ) ‘ INPUT 
READ ( 5 , * ) EF 
WRITE ( 6 , * ) ‘ INPUT 
READ ( 5 , * ) EM 
WRITE (6, *) 1 INPUT 
READ ( 5 , * ) XNUM 
WRITE ( 6 , * ) * INPUT 
READ { 5 , * ) XNUF 
WRITE ( 6 , * ) ' INPUT 
READ ( 5 , * ) RHO 


THE 

THE 

THE 

THE 

THE 

THE 

THE 


COMPOSITE "MATRIX" VOLUME FRACTION' 
COMPOSITE “FIBER" VOLUME FRACTION' 
“FIBER" MODULUS OF ELASTICITY' 
“MATRIX" RELAXATION MODULUS AT t=0 ' 
"MATRIX" POISSONS RATIO AT t=0 ' 
“FIBER" POISSONS RATIO' 

“PLY" DENSITY' 


ENTER NUMBER OF PLIES 


WRITE ( 6 , * ) ' INPUT THE TOTAL NUMBER OF PLIES ' 
READ ( 5 , * ) NLAY 


ENTER THE WIDTH AND LENGTH OF THE BEAM 


WRITE (6, * ) ' INPUT THE WIDTH OF THE CANTILEVER BEAM' 
READ ( 5 , *) WIDTH 

WRITE ( 6 , * ) ' INPUT THE LENGTH OF THE CANTILEVER BEAM' 
READ ( 5 , * ) XLB 
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C 

C CALCULATE SHEAR MODULUS FOR FIBERS AND MATRIX 

C 

GF = EF/ (2* (1+XNUF) ) 

GM = EM/ ( 2 * ( 1+XNUM) ) 


CALCULATE PLY UNIDIRECTIONAL MATERIAL PROPERTIES 
(USING MICROMECHANICS APPROACH) 


DO 10 I = l.NLAY 

E ( 1 , 11) = EF*VF + EM*VM 

E { 1 , 22 ) = EF*EM/(VM*EF + VF*EM) 

GNUM = GF* ( PI+4*VF) + GM*{PI-4*VF) 
GDEN = GF* ( PI-4 *VF) + GM*{PI + 4*VF) 

GN = GNUM /GDEN 

G121 = ( (4-PI)+PI*GN)/4 

G122 = 4*GN/ ( (4-PI) *GN+PI) 

G ( 1 , 12 ) = (GM/2 ) * (G12 1 + G122 ) 

XNU ( I , 12) = XNUF*VF + XNUM*VM 
XNU (1,21) = XNU (I, 12)*E(I,22)/E(I, 11) 
10 CONTINUE 


INPUT THE ANGULAR ORIENTATION OF EACH LAMINA WITH RESPECT TO 
THE LONGITUDINAL AXIS OF THE LAMINATE 


DO 702 IJ = 1 , NLAY 

WRITE ( 6 , * ) 1 INPUT' THE ORIENTATION (ANGLE) OF PLY NO.'.IJ 
READ ( 5 , *) ALPHA ( I J) 

702 CONTINUE 


C 


C 

C 

C 


C 

C 

C 


C 


CALCULATE THE REDUCED STIFFNESS COEFFICIENTS FOR EACH LAMINA 


DO 11 I = 1 , NLAY 

Qll(I) = E ( I , 1 1 ) / ( 1 . -XNU ( I , 12 ) *XNU (1,21)) 

012(1) = XNU (I, 12) *E(I,22)/{1. -XNU ( I , 12 ) *XNU (1,21) ) 
022(1) = E ( I , 22 ) / ( 1 . -XNU ( I , 12 ) *XNU (1,21) ) 

066(1) = G ( I , 12 ) 


CALCULATE THE COMPLIANCE MATRIX COEFFICIENTS FOR EACH LAMINA 


Sll(I) = 1/E( I , 11 ) 

S12 ( I ) = -XNU (I,12)/E(I, 11) 
S22 ( I ) = 1 /E ( I , 22 ) 

S66 ( I ) = 1/G(I, 12) 


CALCULATE THE TRANSFORMED REDUCED STIFFNESSES FOR EACH LAMINA 


ALPHA ( I ) = ALPHA ( I ) *PI /180 . 

ORH(I) = Qll ( I ) * (COS (ALPHA (I) ) ) * *4 

1 +2 . * (Q12 (I)+2 . *Q66 (I) ) * ( (SIN (ALPHA ( I ) ) ) **2) * 

2 ( ( COS ( ALPHA ( I ) ) ) * *2 ) +Q22 ( I ) * (SIN (ALPHA ( I ) ) ) **4 


QR12 ( I ) = (Qll { I ) +Q22 ( I ) -4 . *Q66 ( I ) ) * ( (SIN (ALPHA ( I ) ) ) **2 ) * 
1 ( (COS (ALPHA ( I ) ) ) **2)+Q12 (I) *( (SIN (ALPHA (I ) ) ) **4+ 
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c 


c 


c 


c 


c 

c 

c 


c 


c 


c 


c 


c 


c 

c 

c 

c 

c 

c 


2 ( COS ( ALPHA ( I ) ) ) * * 4 ) 

QRl 6 { I ) = (Qll (I) -Q12 (I) -2 . *Q66 { I ) ) * ( (COS (ALPHA ( I) ) ) **3 ) * 

1 SIN (ALPHA (I))+(Q12(I) -Q22 { I) +2 . *Q66 ( I) ) *COS (ALPHA ( I ) ) * 

2 (SIN (ALPHA ( I ) ) ) **3 


QR2 2(1) = Q11(I)*(SIN(ALPHA(I)))**4 

1 +2.*(Q12(I)+2.*Q66(I) )*( (SIN (ALPHA ( I ) ) ) **2 ) * 

2 ( (COS (ALPHA ( I ) ) ) **2) +Q22 (I) * (COS (ALPHA ( I) ) ) **4 

QR26 ( I ) = (Qll (I) -Q12 (I) -2 . *Q66 ( I) ) * ( (SIN (ALPHA ( I ) ) ) **3 ) * 

1 COS (ALPHA (I))+(Q12(I) -Q22 ( I ) +2 . *Q66 ( I ) ) *S IN (ALPHA ( I ) ) * 

2 (COS (ALPHA ( I ) ) ) **3 

QR66 (I) = (Qll (I)+Q22 (I) -2 . *Q12 ( I) -2 . *Q66 (I) ) * { (SIN (ALPHA ( I ) ) ) **2) 

1 * ( (COS (ALPHA ( I } } ) **2) +Q66 (I) * ( (SIN(ALPHA(I) ) ) **4+ 

2 ( COS ( ALPHA ( I ) ) } * * 4 ) 


CALCULATE THE TRANSFORMED COMPLIANCE COEFFICIENTS FOR EACH LAMINA 


SRll(I) = S 1 1 ( I ) * ( COS ( ALPHA ( I ) ) ) * * 4 

1 + (2 . *S12 ( I ) + S66 ( I) ) * { (SIN (ALPHA ( I ) ) ) **2 ) * 

2 ( (COS (ALPHA (I) ) ) **2) +S22 (I) * (SIN(ALPHA(I) ) ) **4 

SR12 ( I ) = (Sll (I)+S22 (I) -S66 (I) ) * ( (SIN(ALPHA( I ) ) ) **2 ) * 

1 ( (COS ( ALPHA ( I) ) ) **2)+S12 (I) * ( (SIN (ALPHA (I) ) ) **4+ 

2 (COS (ALPHA (I) ) ) **4} 

SRI 6 ( I ) = (2 . *S11 (I) -2 .*S12 (I) -S66 ( I ) ) *( (COS (ALPHA ( I ) ) ) **3 ) * 

1 SIN (ALPHA ( I ) ) - ( 2*S22 (I) -2*S12 (I) -S66 ( I) ) *COS (ALPHA ( I ) 

2 ) * (SIN (ALPHA ( I ) ) ) **3 


SR22 ( I ) = S11(I)*(SIN( ALPHA ( I ) ) ) * * 4 

1 + (2*S12 (I) +S66 (I) ) * ( (SIN (ALPHA ( I ) ) ) **2) * 

2 ( (COS (ALPHA (I) ) )**2)+S22(I)*(COS(ALPHA(I)))**4 

SR26 ( I ) = (2*Sll(I)-2*Sl2(I)-S66(I) )*( (SIN(ALPHA{ I ) ) ) **3 ) * 

1 COS { ALPHA ( I) > - (2*S22 (I) -2*S12 (I) -S66 (I) ) *SIN (ALPHA ( I ) 

2 ) * (COS (ALPHA { I ) ) ) * *3 

SR66 (I) = 2*(2*Sll(I)+2*S22(I)-4 . *S12 (I) -S66 ( I ) ) * ( (SIN (ALPHA ( I ) ) 

1 )**!)*( (COS (ALPHA (I) ) ) **2 ) +S66 ( I ) *( (SIN (ALPHA ( I) ) ) **4+ 

2 ( COS ( ALPHA { I ) ) ) * * 4 ) 

11 CONTINUE 


RELATE EACH LAMINA THICKNESS TO ITS LOCATION WITHIN THE 
LAMINATE (Z COORDINATE IS 0 AT THE MIDSURFACE) 


INPUT LAMINA THICKNESSES 


301 DO 998 J = 1 , NLAY 

IF ( IFLAG . EQ . 0 ) THEN 

WRITE (6,*)' INPUT THE THICKNESS OF PLY NO.‘,J 
READ ( 5 , * ) THK ( J ) 
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END IF 

998 CONTINUE 

IF ( IFLAG . EQ . 1 ) GO TO 302 

C 

C CALCULATE THE WEIGHT OF THE LAMINATE 



TOTWT = 0.0 
DO 700 NA = 1 , NLAY 

TOTWT = TOTWT + WIDTH*THK (NA) *XLB*RHO 
700 CONTINUE 

C 

302 TOTTH =0.0 

DO 25 J = 1 , NLAY 
TOTTH = TOTTH + THK(J) 

25 CONTINUE 

HAFTH = TOTTH/ 2. 

ZT (NLAY) = HAFTH 
DO 24 J = 1 , NLAY 
LAM = NLAY - J + 1 
ZT(NLAY-J) = ZT (LAM) - THK(LAM) 

24 CONTINUE 

IF ( IFLAG . EQ . 1 ) GO TO 300 



C CALCULATE BENDING STIFFNESSES FOR THE ENTIRE LAMINATE 



Dll = 0.0 

D12 =0.0 

D16 =0.0 

D22 =0.0 

D26 =0.0 

D66 =0.0 

DO 35 J = 1 , NLAY 

L = J + 1 

Dll = Dll + QR11 { J) * { (ZT (L-l ) **3 . ) - (ZT (L-2 ) **3 . ) ) /3 . 

D12 = D12 + QR12 ( J) * ( (ZT (L-l ) **3 . ) - ( ZT (L-2 ) **3 . ) ) /3 . 

D16 = D16 + QR16 (J)*((ZT(L-l)**3.)-(ZT(L-2)**3.))/3. 

D22 = D22 + QR22 (J) * ( (ZT(L-l) **3 . ) — ( ZT (L— 2) **3 . ) )/3. 

D26 = D26 + QR26 (J)*( (ZT(L-l) **3 . ) - (ZT (L-2) **3 . ) )/3. 

D66 = D66 + QR66 (J) * { (ZT(L-l) **3 . ) - (ZT(L-2) **3 . ) ) /3 . 

35 CONTINUE 



C CALCULATE THE DETERMINANT OF THE STIFFNESS MATRIX 



DETER = D1 1 *D22 *D66 + D12*D26*D16 + D16*D12*D26 
* - D16*D16*D22 - D12*D12*D66 - D26*D26*Dll 



C CALCULATE COEFFICIENTS FOR THE COFACTORS OF THE STIFFNESS MATRIX 

C 

D11C = D22 *D66 - D26*D26 
D12C = - (D12*D66 - D16*D26) 

D16C = D12 *D26 - D16*D22 
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c 

c 

c 

c 


c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 
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CALCULATE THE COEFFICIENTS OF THE INVERSE STIFFNESS MATRIX 


D11S = D11C/DETER 
D12S = D12C/DETER 
D16S = D16C/DETER 


IDENTIFY APPLIED MOMENTS (APPLIED LOAD) 
(LOAD APPLIED AT TIP [free end] OF BEAM) 


WRITE (6,*) ' INPUT APPLIED LOAD AT THE TIP OF THE CANTILEVER BEAM 1 

READ ( 5 , * ) PLOAD 

XMAP = - PLOAD *XLB /WIDTH 


*************** PERFORM DEFLECTION CALCULATIONS ****************** 


CALCULATE THE CURVATURE AT ANY POINT ALONG THE LENGTH OF THE BE AM 


INMl = 30 
XIN(l) = 0.0 
DO 805 K = 2,31 
XIN(K) = XIN(K-l) +XLB/INM1 
805 CONTINUE 

DO 804 NI = 1 , INM1+1 

CAL(NI) = DllS*PLOAD* (XLB-XIN(NI) ) /WIDTH 
804 CONTINUE 


PERFORM A POLYNOMIAL CURVE FIT FOR THE CURVATURE AT 
ANY POINT ALONG THE LENGTH 


CALL POLYCF (XIN, CAL, 31,3, XL AMD, NERR) 


CALCULATE THE TIP DEFLECTION OF THE BEAM BY INTEGRATING THE 
CURVATURE TWICE WITH RESPECT TO THE LENGTH COORDINATE 


DEF =0.0 
DO 111 J = 1, 4 

DEF = DEF + XLAMD { J) * (XLB* * { J+l ) ) / (J* (J+l) ) 

111 CONTINUE 

WRITE (13, 109) DEF 

109 FORMAT ( 1 THE TIP DEFLECTION AT TIME ZERO IS , ,E12.5) 


CALCULATE THE NATURAL FREQUENCY OF THE CANTILEVER BEAM 


FREQ = 3 . 51601 *SQRT ( (WIDTH *3 86 . 4 ) / (DllS*TOTWT* (XLB**4 ) ) ) 


WRITE { 13 , * ) * D11S =■ , DllS 
WRITE { 13 , * ) 1 WIDTH = 1 , WIDTH 
WRITE { 13 , * ) 1 WEIGHT = 1 , TOTWT 
WRITE ( 13 , * ) ' LENGTH = ' , XLB 
WRITE (13, 701) FREQ 

7 01 FORMAT (' THE NATURAL FREQUENCY AT TIME ZERO IS',E12.5, ' HERTZ') 
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CALCULATE THE STRESSES FOR EACH LAMINA 


I FLAG = 1 
GO TO 301 

300 DO 36 J = 1 , NLAY 
L = J + 1 - 2 

SIGX(J) = ZT (L) *XMAP* (QR1 1 ( J) *D1 IS +QR12 (J) *D12S +QR16 (J) *D16S) 
SIGY(J) = ZT (L) *XMAP* (QR12 ( J) *D11S +QR22 ( J) *D12S +QR26 { J) *D16S) 
SIGXY(J) = ZT (L) *XMAP* (QR16 { J) *D11S +QR26 ( J) *D12S +QR66 ( J) *D16S) 

36 CONTINUE 

C 

C CALCULATE THE INITIAL (T=0.0) STRAINS FOR EACH LAMINA 

C 

WRITEQ3, 251) 

251 FORMAT ( 1 PLY ' , 7X, 1 EPSX ' , 11X, 1 EPSY MIX, • EPSXY ' / ) 

DO 37 J = 1 , NLAY 

EPSX ( J) = SRI 1 ( J ) *SIGX ( J ) +SR12 ( J) *SIGY { J) +SR16 ( J) *SIGXY ( J) 
EPSY(J) = SR 12 ( J) *SIGX ( J) +SR22 ( J ) *SIGY ( J) +SR26 { J) *SIGXY ( J) 

EPSXY ( J ) = SR 16 (J) *SIGX (J) +SR26 { J) *SIGY ( J) +SR66 (J) *SIGXY (J) 

WRITE (13, 250) J, EPSX (J) , EPSY(J) , EPSXY ( J) 

250 FORMAT ( 1 4 , 3 E 1 5 . 6 ) 

37 CONTINUE 

C 

STOP 

END 
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